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ABSTRACT 


Many  flexible  multibody  applications  are  characterized  by  high  inertia  forces  and  motion 
discontinuities.  Because  of  these  characteristics,  problems  can  be  encountered  when  large 
displacement  finite  element  formulations  are  used  in  the  simulation  of  flexible  multibody 
systems.  In  this  investigation,  the  performance  of  two  different  large  displacement  finite 
element  formulations  in  the  analysis  of  flexible  multibody  systems  is  investigated.  These  are 
the  incremental  corotational  procedure  proposed  by  Rankin  and  Brogan  [15]  and  the  non- 
incremental  absolute  nodal  coordinate  formulation  recently  proposed  [19].  It  is  demonstrated 
in  this  investigation  that  the  limitation  resulting  from  the  use  of  the  nodal  rotations  in  the 
incremental  corotational  procedure  can  lead  to  simulation  problems  even  when  very  simple 
flexible  multibody  applications  are  considered.  The  absolute  nodal  coordinate  formulation, 
on  the  other  hand,  does  not  employ  infinitesimal  or  finite  rotation  coordinates  and  leads  to 
a  constant  mass  matrix.  Despite  the  fact  that  the  absolute  nodal  coordinate  formulation 
leads  to  a  complex  expression  for  the  elastic  forces,  the  results  presented  in  this  study, 
surprisingly,  demonstrate  that  such  a  formulation  is  efficient  in  static  problems  as  compared 
to  the  incremental  corotational  procedure.  The  excellent  performance  of  the  absolute  nodal 
coordinate  formulation  in  static  and  dynamic  problems  can  be  attributed  to  the  fact  that 
such  a  formulation  does  not  employ  rotations  and  leads  to  exact  representation  of  the  rigid 
body  motion  of  the  finite  element. 


1  INTRODUCTION 


The  computational  issues  associated  with  large  displacement  problems  [2,  4,  6,  10,  16,  20] 

.  become  important  when  flexible  multibody  applications  are  considered.  This  is  due  to  the 
high  nonlinearities  in  the  equations  of  motion,  the  coupling  of  the  elastic  and  reference 
motions,  high  inertia  forces  and  possible  motion  discontinuities.  Therefore,  it  is  important 
to  carefully  examine  the  accuracy,  robustness  and  efficiency  of  the  computational  procedures 
used  in  the  large  displacements  of  multibody  system  applications.  The  most  widely  used 
formulation  in  flexible  multibody  dynamics  is  the  floating  frame  of  reference  formulation 
[12,  13,  14,  19].  The  use  of  this  formulation,  however,  has  been  limited  to  small  deformation 
problems.  In  the  floating  frame  of  reference  formulation,  the  body  elastic  deformation  is 
described  in  a  body  coordinate  system,  and  it  is  assumed  that  this  deformation  is  small 
e  in  order  to  justify  the  use  of  linear  modes  [19].  Consequently,  this  formulation  has  been 

rarely  used  for  large  deformation  problems,  which  can  be  analyzed  more  accurately  using  a 

* 

full  finite  element  representation.  Nonetheless,  large  deformation  problems  can  be  examined 
using  the  floating  frame  of  reference  formulation  by  dividing  the  flexible  body  into  a  large 
number  of  bodies,  each  of  which  has  its  own  body  reference.  This  approach,  however,  leads 
to  large  dimensionality  and  nonlinearity  of  the  inertia  forces. 

There  are  two  types  of  finite  element  procedures  that  can  be  used  for  the  large  defor¬ 
mation  analysis;  incremental  and  non-incremental.  The  incremental  approach  is  the  most 
widely  used  procedure  for  the  solution  of  non-linear  large  rotations  and  large  deformation 
problems  in  structural  applications.  Several  incremental  procedures  have  been  developed  for 
non-isoparametric  elements,  in  which  infinitesimal  rotations  are  used  as  nodal  coordinates. 
In  principle,  these  procedures  can  also  be  used  in  multibody  applications,  since  such  proce¬ 
dures  can  also  be  used  to  describe  the  large  reference  displacements  and  rotations  which  are 
s  characteristics  of  multibody  problems.  However,  the  limitations  on  the  rotation  increments 

in  some  of  these  procedures,  as  will  be  discussed  in  this  paper,  and  the  fact  that  these  pro- 
i 
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cedures  do  not  lead  to  an  exact  rigid  body  inertia  representation  as  the  result  of  the  early 
linearization  of  the  equations  of  motion  [17],  make  the  incremental  approach  less  attractive 
to  use  in  .flexible  multibody  problems.  Belytschko  and  Hsieh  [6]  used  the  convected  coordinate 
system  and  applied  it  to  the  dynamic  analysis  of  structural  systems  that  undergo  large  rota¬ 
tions.  In  this  formulation,  a  convected  coordinate  system  is  assigned  to  each  finite  element 
and  the  element  internal  forces  are  first  defined  in  the  element  convected  coordinate  system 
and  then  transformed  to  the  global  coordinate  system.  A  basic  feature  of  this  technique 
is  the  decomposition  of  the  global  displacement  field  into  rigid-body  and  strain-producing 
deformation  components.  The  increment  steps  are  chosen  such  that  the  element  rotation 
between  two  consecutive  configurations  is  small  and  the  element  shape  function  and  local 
nodal  coordinates  can  be  used  to  describe  this  small  rotation.  Argyris  et  al.  [2]  presented  a 
detailed  discussion  on  the  convected  coordinate  procedure  and  the  large  deflection  problems. 
They  introduced  the  natural  approach  that  refers  to  the  separation  between  the  rigid  body 
displacement  field  and  the  natural  deformation  in  the  total  displacement  field  of  a  finite 
element.  Hughes  and  Winget  [10]  presented  an  efficient  algorithm  to  define  the  displace¬ 
ment  increment  over  the  step  in  the  large  deformation  analysis  and  demonstrated  that  a 
unique  large  rotation  vector  can  be  assigned  to  any  rotation  in  large  deformation  problems. 
A  corotational  procedure  for  the  solution  of  nonlinear  finite  element  large  rotation  problems 
was  proposed  by  Rankin  and  Brogan  [15].  This  procedure  will  be  discussed  in  detail  in  the 
following  section  and  will  be  used  in  the  study  presented  in  this  paper.  In  this  procedure, 
the  contribution  of  the  so  called  large  rigid-body  rotations  of  the  element  is  removed  from 
the  global  displacement  field  through  the  use  of  an  element  convected  coordinate  system. 
A  nonsingular  large  rotation  vector  is  introduced  to  describe  the  nodal  rotations  [3],  Hsiao 
and  Jang  [9]  extended  the  use  of  the  corotational  procedure  to  the  dynamic  analysis  of  pla¬ 
nar  flexible  linkages.  A  detailed  corotational  formulation  for  the  dynamic  analysis  of  planar 
beams  undergoing  large  deflections  has  been  recently  presented  by  Behdinan  et  al.  [5] . 
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A  new  non-incremental  approach,  the  absolute  nodal  coordinate  formulation,  has  recently- 
been  proposed  [19].  This  formulation  differs  from  other  existing  finite  element  formulations 
in  the  sense  that  no  infinitesimal  or  finite  rotations  are  used  as  nodal  coordinates.  The 
set  of  nodal  coordinates  consists  of  global  displacements  and  slopes.  Using  this  approach, 
beams  and  plates  can  be  treated  as  isoparametric  elements,  and  therefore  there  is  no  need 
to  introduce  an  element  coordinate  system  to  describe  the  rigid  body  rotations  of  the  finite 
element.  The  absolute  nodal  coordinate  formulation  leads  to  a  constant  mass  matrix,  while 
the  elastic  forces  are  nonlinear  functions  of  the  element  coordinates.  In  the  absolute  nodal 
coordinate  formulation,  large  rigid  body  displacements  including  large  rotations  produce  zero 
strains  in  the  finite  elements.  It  was  demonstrated  [19]  that  in  order  to  obtain  correct  results 
in  the  dynamic  analysis,  a  consistent  mass  approach  must  be  used  in  this  formulation. 

It  is  the  objective  of  this  paper  to  examine  the  performance  of  the  absolute  nodal  co¬ 
ordinate  formulation  by  comparing  it  with  the  corotational  procedure  presented  by  Rankin 
and  Brogan  [15]  and  implemented  in  the  finite  element  code  ANSYS  [1].  It  will  be  shown 
that  numerical  problems  are  encountered  when  the  incremental  procedure  is  used  in  flexible 
multibody  applications.  This  paper  is  organized  as  follows.  In  Section  2  the  incremental 
corotational  procedure  proposed  by  Rankin  and  Brogan  [15]  is  reviewed.  This  procedure  will 
be  extensively  used  in  our  investigation,  and  therefore,  the  review  materials  presented  in  Sec¬ 
tion  2  are  used  as  the  basis  for  the  discussion  presented  in  the  following  sections.  In  Section 
3,  the  non-incremental  absolute  nodal  coordinate  formulation  is  introduced.  In  Section  4,  the 
performance  of  the  absolute  nodal  coordinate  formulation  in  static  problems  is  examined  and 
compared  with  the  incremental  procedure.  Two  problems  with  known  analytical  solutions 
are  considered.  These  are  the  elastica  problem  and  the  bending  of  a  beam  into  a  full  circle. 
In  Section  5,  the  main  features  of  the  absolute  nodal  coordinate  formulation  in  the  case  of 
dynamics  are  summarized.  Comparison  between  the  incremental  corotational  procedure  and 
the- non-incremental  absolute  nodal  coordinate  formulation  when  flexible  multibody  appli- 


cations  are  considered  is  presented  in  Section  6.  Summary  and  conclusions  drawn  from  this 
study  are  presented  in  Section  7. 


2  COROTATIONAL  PROCEDURE 

The  incremental  procedure  has  been  widely  and  successfully  used  in  the  nonlinear  finite 
element  analysis  of  large  rotation  structural  problems.  The  incremental  finite  element  coro- 
tational  procedure  proposed  by  Rankin  and  Brogan  [15]  has  been  implemented  in  several 
general  purpose  structural  analysis  codes  such  as  ANSYS,  and  has  been  used  in  the  analysis 
of  many  large  rotation  and  deformation  problems.  In  this  procedure,  which  is  independent 
of  the  element  formulation,  any  rigid  body  motion  contribution  is  eliminated  from  the  global 
displacement  field  in  order  to  determine  the  pure  deformation.  The  contribution  of  the  rigid 
body  rotations  of  the  element  is  eliminated  by  using  a  convected  coordinate  system  that 
moves  with  the  element.  The  element  equations  are  first  defined  in  the  element  coordinate 
system  and  then  transformed  in  order  to  define  these  equations  in  the  global  inertial  frame 
These  equations  are  solved  for  the  displacement  increments  that  are  then  used  to  update  the 
global  displacement  field  of  the  element. 

In  this  approach,  the  nonlinear  kinematics  of  the  finite  element  is  defined  in  terms  of  a 
large  reference  motion  plus  a  small  deformation;  this  holds  assuming  that  at  each  time  step 
the  displacement  increments  are  so  small  that  the  current  configuration  in  which  the  element 
equations  are  defined  can  be  considered  a  valid  reference  configuration.  This  implies  that  in 
one  time  step  there  is  no  large  variation  in:  1)  the  deformation  within  each  element,  and  2) 
the  large  reference  motion.  Consequently,  the  most  important  parameter  that  governs  this 
procedure  is  the  time/load  step  which  must  remain  small.  It  will  be  shown  later  that  there 
is  another  limitation  due  to  the  assumption  that  the  total  deformation  within  each  element 
must  remain  small. 


* 

In  order  to  extract  the  rigid  body  motion,  a  local  coordinate  system  is  introduced.  Using 
the  notation  of  Rankin  and  Brogan  [15],  let  Ek  be  the  orthogonal  transformation  matrix 
that  defines  the  orientation  of  the  local  element  frame  in  the  global  frame  at  the  k- th  step. 
A  rigid  body  rotation  can  be  extracted  from  the  total  displacement  as  follows.  The  portion 
ufe^  of  the  total  displacement  that  causes  strain  is  given  by 

uf'  =  Ej(u,  +  Xs)-Xe,  (1) 

where  us  is  the  total  displacement  defined  in  the  inertial  coordinate  system,  X5  is  the  global 
position  of  an  arbitrary  point,  and  Xe  is  the  local  position  of  the  same  point  before  defor¬ 
mation.  The  vector  ufCjf  will  be  used  to  define  the  strain  energy  and  the  generalized  elastic 
forces.  The  solution  of  the  system  of  equations  at  step  k  yields  a  displacement  increment 
Aufc+i.  Using  this  increment,  it  is  possible  to  calculate  the  total  displacement  ufc+1  and  use 

m) 

it  to  define  the  new  orthogonal  transformation  matrix  Efc+1.  The  small  rotation  increment  is 
«.  used  to  update  the  rotations  within  the  element  using  the  corotational  approach.  In  fact,  the 

deformational  rotations  can  be  finite  rotations,  and  cannot  be  treated  as  ordinary  vectors 
[3].  For  this  reason,  the  rotations  are  described  in  terms  of  pseudovectors  and  the  nodal 
deformation  rotational  degrees  of  freedom  are  treated  differently  from  nodal  deformation 
translational  degrees  of  freedom.  Rankin  and  Brogan  [15]  introduced  a  surface  coordinate 
system  rigidly  attached  to  each  node.  This  surface  coordinate  system  is  defined  in  the  inertial 
frame  by  the  orthogonal  transformation  Sk,  as  shown  in  Fig.  1.  Clearly  the  deformation  is 
produced  by  a  relative  rotation  of  Sjt  with  respect  to  E;. .  This  relative  rotation  is  expressed 
through  the  orthogonal  matrix  T^,  where 

T*  =  Ej'StSj'Eo.  (2) 

Hughes  and  Winget  [10]  have  shown  that  the  quantity 

* 

fl  —  2(T  —  I)(T  + 1)-1  (3) 
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is  always  skew-symmetric  for  any  orthogonal  matrix  T.  Rankin  and  Brogan  [15]  first  consid¬ 
ered  the  following  definition  of  ft,  that  is  slightly  different  from  the  one  presented  by  Argyris 
[3]: 

0  -u3  U>2 

ft  =  to3  0  —to  i 

—U>2  tO\  0 

where  7  is  the  trace  of  T: 

7  —  tr(T)  =  1  +  2cos  9,  (5) 


2(T  -  TT) 
1  +  7  ’ 


(4) 


and  6  is  the  angle  of  rotation  about  the  principal  eigenvector  of  T.  The  matrix  ft  is  associated 
with  the  vector 


to  = 


T 


LO\  CJ2  ^3 


(6) 


which  is  called  rotational  pseudovector ,  and  becomes  the  rotation  vector  when  rotations  are 
small.  Given  the  pseudovector  to,  the  skew-symmetric  matrix  is  calculated  and  used  to 
evaluate  the  matrix  T  as 


T  =  I- 


'i  +  i 


\to\ 


Consequently,  the  norm  |u?|  of  to  is  related  to  6  as 


(7) 


|u>|  =  2  tan 


(8) 


In  order  to  avoid  singularity,  Rankin  and  Brogan  use  a  different  definition  of  the  transfor¬ 
mation  matrix,  that  is  given  by 

T  =  I  +  fi^l-iM!+^2,  (9) 

and,  in  this  case,  the  relationship  between  |u;|  and  6  is 

M  =2sin^.  (10) 
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This  equation  shows  the  limit  of  Rankin  and  Brogan’s  corotational  procedure.  Using  the 
elements  of  oj  as  rotational  degrees  of  freedom,  the  magnitude  |u>|  can  approximate  the  angle 
9  only  for  rotations  up  to  30°. 

The  algorithm  proposed  by  Rankin  and  Brogan  can  be  summarized  as  follows.  Given  the 
configuration  of  a  finite  element  at  step  k,  local  translational  displacements  are  calculated 
using  Eq.  1.  A  relative  rotation  matrix  Tfc  is  calculated  using  Eq.  2.  The  rotational  freedoms 
are  taken  from  the  pseudovector  that  is  associated  with  the  skew-symmetric  matrix 


^ k 


2(T»  -  Tf) 


(11) 


1  +  tr(Tfc) 

Internal  forces  calculated  within  the  element  are  then  transformed  to  the  global  coordinate 
system.  Displacement  increments  required  to  reduce  the  out-of-balance  forces  are  then  ob¬ 
tained  by  solution  of  the  system  of  equations.  While  the  translational  increments  are  directly 
added  to  the  total  displacement  vector,  the  rotational  increments  Au:  are  used  to  update 
the  nodal  pseudovector  u>^  using  the  following  equation: 

Auj^Jl  -  \  |a?fc|2  +  LQk  ~  x  Aw 

\A  +  \  |o>fc|2 

where  the  sign  must  be  the  same  as  the  sign  of  the  quantity 


Wfc+i  =  ±- 


(12) 
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'1  -  j  Kf Au>. 


(13) 


Note  that  the  pseudovector  cok  defines  the  relative  rotation  of  the  surface  coordinate  system 
with  respect  to  the  element  coordinate  system.  The  global  rotations  are  not  updated  using 
the  pseudovector  approach,  but  the  rotation  increments  are  simply  added  to  the  rotations 
at  the  previous  step.  A  complete  discussion  on  the  corotational  procedure  can  be  found  in 
[1]  and  [15]. 

As'  previously  pointed  out,  the  corotational  procedure  has  been  widely  used  and  imple¬ 
mented  in  general  purpose  finite  element  codes  (e.g.  STAGS,  ANSYS).  It  was  demonstrated 
that  this  procedure  is  efficient  and  accurate  in  large  rotation/small  strain  problems.  When 


the  strain  within  the  element  becomes  larger,  this  procedure  does  not  perform  well.  It  was 
shown  in  this  section  that  the  element  local  rotations  (the  rotations  of  the  surface  coordinate 
system  relative  to  the  element  coordinate  system)  must  remain  less  than  30°,  otherwise  the 
approximation  of  the  element  rotational  degrees  of  freedom  with  the  pseudovector  compo¬ 
nents  can  not  be  considered  accurate.  Furthermore,  in  order  to  define  correctly  the  element 
equations  in  the  current  element  configuration,  the  element  displacement  increments  must 
be  small.  While  these  restrictions  may  not  be  serious  in  structural  applications,  especially 
when  very  small  load/ time  steps  are  used,  serious  problems  can  be  encountered  when  simu¬ 
lating  simple  dynamic  multibody  problems  as  demonstrated  by  the  results  presented  in  later 
sections. 


3  NON  INCREMENTAL  ABSOLUTE  NODAL  COORDINATE  FORMULA¬ 
TION 

In  this  section,  the  non-incremental  finite  element  absolute  nodal  coordinate  formulation  [19] 
is  briefly  reviewed.  In  this  formulation,  the  nodal  coordinates  of  the  element  are  defined  in 
a  fixed  inertial  coordinate  system,  and  consequently  no  transformation  is  required  for  the 
element  coordinates.  The  element  nodal  coordinates  represent  global  nodal  displacements 
and  slopes.  Thus,  in  the  absolute  nodal  coordinate  formulation,  no  infinitesimal  or  finite 
rotations  are  used  as  nodal  coordinates  and  no  assumption  on  the  magnitude  of  the  element 
rotations  is  made.  The  results  presented  in  this  paper  will  demonstrate  that  these  properties 
of  the  absolute  nodal  coordinate  formulation  make  this  formulation  efficient  and  accurate  in 
the  large  displacement  analysis  of  flexible  multibody  systems. 

In  this  investigation,  two  dimensional  beam  elements  are  considered.  The  global  position 
vector  r  of  an  arbitrary  point  P  on  the  element  is  defined  in  terms  of  the  nodal  coordinates 
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and  the  element  shape  fimction,  as  shown  in  Fig.  2,  as 


r  = 


n 

r2 


=  Se, 


(14) 


where  S  is  the  global  shape  function  which  has  a  complete  set  of  rigid  body  modes,  and  e  is 
the  vector  of  element  nodal  coordinates: 


e  = 


lT 


ei  e2  63  e4  es  e-j  eg 


(15) 


This  vector  of  absolute  nodal  coordinates  includes  the  global  displacements 


ei  =  r!|x=0 ,  e2  =  r2\x=0 ,  e5  =  rx\x=l ,  e6  =  r2|x==/ , 
and  the  global  slopes  of  the  element  nodes,  that  are  defined  as 


(16) 


drx 

dr2 

drx 

dr2 

dx 

3 

a;=0 

64 

3 

x—0 

>  =  -W~ 

_  ,  OX 

x—L 

(17) 


x=l 


Here  x  is  the  coordinate  of  an  arbitrary  point  on  the  element  in  the  undeformed  con  figuration 
and  l  is  the  original  length  of  the  beam  element.  Since  absolute  coordinates  are  used,  a  cubic 
polynomial  is  employed  to  describe  both  components  of  the  displacements.  Therefore,  the 
global  shape  function  S  can  be  written  as 


si  0  s2l  0  S3  0  s^l  0 

0  Si  0  s2l  0  s3  0  S4I 

where  the  functions  Sj  =  s*(£)  are  defined  as 


(18) 


si  =  1  —  3£2  +  2£3,  s2=e- 2£2  +  £3,  s3  =  3f2-2f3,  s4  =  f3-^,  (19) 

and  £  =  x/l.  It  can  be  shown  that  the  preceding  shape  function  contains  a  complete  set 
of  rigid  body  modes  that  can  describe  arbitrary  rigid  body  translational  and  rotational 
displacements,  provided  that  global  slope  coordinates  are  used  instead  of  infinitesimal  rota¬ 
tions.  Using  the  absolute  coordinates  and  slopes,  it  can  also  be  shown  that  the  beam  element 
defined  by  the  shape  function  of  Eq.  18  is  an  isoparametric  element. 
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4  PERFORMANCE  IN  STATIC  PROBLEMS 

In  this  section,  the  performance  of  the  two  incremental  and  non-incremental  finite  element 
procedures  discussed  in  the  preceding  two  sections  in  the  static  analysis  of  large  deflection 
problems  of  planar  beams  is  investigated  using  two  numerical  examples.  These  two  examples 
of  large  static  deformation  are  solved  using  the  general  purpose  finite  element  code  ANSYS 
that  utilizes  the  corotational  procedure  presented  by  Rankin  and  Brogan  [15]  for  the  large 
deformation  problems  [1],  and  they  are  also  solved  using  the  non-incremental  absolute  nodal 
coordinate  formulation  [7].  The  results  of  ANSYS  and  the  absolute  nodal  coordinate  formu¬ 
lation  are  obtained  using  a  linear  strain-displacement  relationship.  The  first  example  is  a 
cantilever  beam  loaded  with  a  free  end  moment  that  bends  into  a  full  circle,  while  the  second 
problem  is  the  elastica  problem.  Both  examples,  which  have  a  known  analytical  solution, 
employ  the  same  beam  model.  The  beam  in  this  model  is  assumed  to  have  length  of  1  m, 
cross  sectional  area  of  1.257E-03  m2,  second  moment  of  area  of  1.257E-07  m4,  and  modulus 
of  elasticity  of  2.0E+09  Pa. 

Bending  of  a  Beam  into  a  Full  Circle  This  example  is  shown  in  Fig.  3.  The  beam 
is  divided  into  10  elements  and  the  results  obtained  using  ANSYS  are  almost  identical  to  the 
exact  solution.  The  total  CPU  time  required  to  obtain  the  solution  shown  in  Fig.  3  using 
ANSYS  on  HP-Convex  SPP1200/XA-16  was  found  to  be  7.5  sec.  The  same  results  were  also 
obtained  by  Rankin  and  Brogan,  who  emphasized  that  the  corotational  procedure  gives  good 
results  in  this  kind  of  analysis  because  the  large  deflections  of  the  beam  are  converted  into 
much  smaller  deformational  increments  at  each  load  step.  They  also  showed  that  the  results 
obtained  using  a  conventional  incremental  approach  were  not  accurate  and  that  no  solution 
could  be  obtained  when  the  free  end  rotation  reached  about  90°  [15].  In  the  example  of 
Fig.  3i  the  ratio  between  the  nodal  rotation  increments  and  the  load  increments  is  constant. 
Furthermore,  the  curvature  of  the  beam  remains  constant  along  its  length  as  shown  in  Fig. 
3.  Because  of  these  characteristics,  the  load  step  increments  are  easily  generated  such  that 
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the  rotational  increments  lie  within  the  range  of  allowed  rotations. 

The  absolute  nodal  coordinate  formulation  was  also  used  to  solve  the  same  problem. 
Figure  4  shows  the  results  of  the  global  rotation  of  the  beam  free  end  in  the  full  circle 
example.  The  results  obtained  using  the  absolute  nodal  coordinate  formulation  are  compared 
with  the  exact  solution.  In  Fig.  5,  the  solution  configurations  of  the  beam  loaded  by  the 
free  end  moment  are  presented.  These  results  were  obtained  by  dividing  the  beam  into  10 
finite  elements.  The  total  CPU  time  for  obtaining  this  solution  on  a  PC  Pentium  90MHz 
was  found  to  be  7.2  sec. 

Figure  4  shows  that  the  absolute  nodal  coordinate  formulation  gives  slightly  different 
results  from  the  exact  solution  (6%  error)  when  9}reeend  =  360°.  Nonetheless,  the  solution 
obtained  using  the  absolute  nodal  coordinate  formulation  is  quite  accurate  and  the  CPU 
time  demonstrates  that  the  method  is  computationally  efficient.  The  difference  from  the 
exact  solution  is  due  to  the  fact  that  a  linear  elastic  model  was  used  for  the  formulation  of 
the  elastic  forces.  It  can  be  demonstrated  that  a  better  solution  is  obtained  using  a  higher 
number  of  elements. 

Elastica  Problem  The  second  example  of  static  analysis  of  beam  large  deflection 
problem  is  the  elastica  problem  in  which  a  cantilever  beam  is  subject  to  a  compressive  load 
at  the  free  end  [7].  The  analytical  solution  of  this  problem  can  be  found  in  Timoshenko 
[21].  Figure  6  shows  the  deformed  configurations  of  the  beam  predicted  using  ANSYS  under 
different  loads  over  the  critical  limit  for  the  cases  of  10  and  20  finite  elements  discretization. 
In  Table  1,  the  results  of  the  global  rotations  of  the  free  end  node  obtained  using  ANSYS  are 
presented  and  compared  to  the  exact  solutions.  The  CPU  time  for  obtaining  the  solution  of 
Fig.  6  on  an  HP-Convex  SPP1200/XA-16  was  60  sec  and  64  sec  for  10  and  20  finite  element 
cases,  respectively.  The  solution  obtained  using  ANSYS  is  very  close  to  the  exact  solution 
in  the  range  40°-120°of  rotation  of  the  free  end.  Rankin  and  Brogan  [15]  affirmed  that 
no  solution  to  the  elastica  problem  could  be  obtained  using  the  conventional  incremental 
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approach.  When  the  load  approaches  the  critical  value  (for  free  end  rotations  <  40°),  the 
results  from  the  ANSYS  solution  are  different  from  the  exact  solution.  This  is  due  to  the 
fact  that  in  the  vicinity  of  the  critical  value,  a  very  small  change  of  load  causes  a  very  large 
change  in  solution  configurations.  In  fact,  when  the  tolerance  parameter  for  the  equilibrium 
iterations  is  decreased  for  a  better  iterative  refinement,  a  solution  much  closer  to  the  exact 
one  can  be  obtained  in  the  range  of  solutions  with  free  end  rotations  <  40°.  However,  it  was 
impossible  to  obtain  a  solution  using  ANSYS  for  loads  that  give  free  end  rotations  >  140° 
using  both  10  and  20  elements.  The  solution  for  a  load  which  causes  an  exact  solution  of 
140°  of  free  end  rotation  is  also  inaccurate.  Close  analysis  of  this  problem  showed  that  for 
loads  that  give  free  end  rotations  >  140°  the  curvature  of  the  beam  in  the  neighborhood 
of  the  fixed  end  becomes  quite  large,  while  the  rest  of  the  beam  remains  almost  straight. 
When  the  curvature  becomes  large,  the  relative  rotation  is  greater  than  30°.  Furthermore, 
the  automatic  load  stepping  is  governed  by  the  average  element  rotational  increments  that 
in  this  example  remain  small  even  though  the  curvature  in  the  neighborhood  of  the  fixed 
end  is  relatively  large.  Consequently,  it  becomes  very  difficult  for  the  ANSYS  code  to  adjust 
the  load  step  according  to  the  increase  of  the  curvature  in  a  very  limited  area  of  the  beam 
This  example  shows  that  using  the  rotations  as  nodal  degrees  of  freedom  leads  eventually  to 
accuracy  and  convergence  problems. 

The  same  elastica  problem  was  solved  using  the  absolute  nodal  coordinate  formulation. 
Table  2  shows  the  results  of  the  global  rotation  of  the  beam  free  end  in  the  elastica  problem 
using  the  absolute  nodal  coordinate  formulation  with  10  and  20  finite  elements  [7].  The 
solution  configurations  under  different  overcritical  loads  are  shown  in  Fig.  7.  These  results 
show  that  the  solution  obtained  using  the  absolute  nodal  coordinate  formulation  is  very 
close  to  the  exact  solution.  For  loads  close  to  the  critical  value  (solutions  with  free  end 
rotations  <  40°),  the  case  of  10-element  discretization  gives  better  results  than  the  20- 
element  case,  because  the  solution  is  obtained  for  a  smaller  number  of  variables  while  the 
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elastic  deformation  of  the  beam  remains  relatively  small.  Unlike  the  results  obtained  using 
ANSYS,'  accurate  solutions  are  also  obtained  for  free  end  rotations  >  140°.  In  these  cases, 
the  20-element,  model  performs  better  than  the  10-element  model,  since  the  more  refined 
discretization  of  the  beam  can  better  describe  the  large  deformations.  The  total  CPU  time 
used  to  obtain  the  solution  on  a  Pentium  90  MHz  was  8  sec  and  14  sec  for  the  10  and  20- 
element  models,  respectively.  These  results  demonstrate  that  the  non-incremental  absolute 
nodal  coordinate  formulation  leads  to  an  accurate  and  efficient  solution  as  compared  to  the 
incremental  corotational  procedure. 


5  DYNAMIC  PROBLEMS 

The  results  presented  in  the  preceding  section  demonstrate  that  the  non-incremental  absolute 
nodal  coordinate  formulation  performs  well  in  static  problems  despite  the  fact  that  such  a 
formulation  leads  to  a  complex  expression  for  the  elastic  forces.  In  fact,  it  is  surprising  to 
note  that  the  absolute  nodal  coordinate  formulation  is  more  efficient  as  compared  to  the 
incremental  methods  in  static  applications,  and  this  formulation  leads  to  accurate  results  in 
the  large  deformation  problems  by  using  a  linear  strain-displacement  relationship.  Since  the 
absolute  nodal  coordinate  formulation  leads  to  a  constant  mass  matrix,  it  is  expected  that 
this  formulation  will  perform  even  better  in  dynamics  problems. 

Incremental  Finite  Element  Approach  The  incremental  finite  element  approach 
has  been  widely  used  for  the  dynamic  analysis  of  flexible  systems  that  undergo  large  rotations 
and  deformations.  In  the  incremental  finite  element  formulation,  nodal  rotations  are  used 
as  degrees  of  freedom  and  the  nodal  coordinates  are  treated  as  vectors  [2,  10,  15].  The 
internal  forces  of  the  flexible  bodies  are  first  defined  in  the  element  coordinate  systems  and 
then  transformed  to  the  global  system.  The  dynamic  equations  are  then  solved  for  the 
deformation  increments.  In  Section  2,  a  corotational  procedure  was  presented  which  allows 
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the  use  of  the  conventional  finite  element  formulations  in  large  rotation  problems.  In  the 
dynamics  of  flexible  bodies  that  undergo  large  rotations,  it  is  important  to  obtain  accurate 
modeling  of  the  inertia  of  the  bodies.  However,  it  was  recently  demonstrated  [17]  that  the 
use  of  the  incremental  approach  where  rotations  are  used  as  nodal  coordinates  does  not  lead 
to  the  exact  modeling  of  the  rigid  body  dynamics  of  simple  structures. 

When  the  incremental  formulations  are  used  with  consistent  mass  techniques,  the  global 
mass  matrix  of  the  element  is  not  constant.  As  a  consequence,  the  expression  of  the  inertia 
forces  does  not  take  a  simple  form  and  these  forces  have  to  be  updated  at  every  time  step. 
In  the  ANSYS  code,  the  conventional  shape  function  of  the  beam  element  is  used.  In  this 
shape  function,  the  axial  displacement  is  approximated  using  a  linear  polynomial,  while  the 
transverse  displacement  is  approximated  using  a  cubic  polynomial.  This  is  the  displacement 
field  which  is  used  to  generate  the  ANSYS  results  presented  in  the  following  section. 

Absolute  Nodal  Coordinate  Formulation  In  Section  3,  the  generalized  nodal  co¬ 
ordinates  and  the  displacement  field  of  the  absolute  nodal  coordinate  formulation  were  pre¬ 
sented.  In  this  formulation,  the  global  position  vector  of  an  arbitrary  point  on  the  element 
is  defined  in  terms  of  a  set  of  global  nodal  coordinates  and  a  global  shape  function.  It  is 
assumed  that  this  global  shape  function  has  a  complete  set  of  rigid  body  modes.  By  differ¬ 
entiating  Eq.  14  with  respect  to  time,  we  obtain  the  global  velocity  vector  that  can  be  used 
to  define  the  kinetic  energy  of  the  element  as 

T=IJV  fiTidV  =  bT  (X  »STSdl/)  e,  (20) 

where  p  and  V  are,  respectively,  the  mass  density  and  volume  of  the  element.  We  can  define 
the  mass  matrix  of  the  element  as 

M  =  JvpSrSdV,  (21) 

where  M  is  a  symmetric  and  constant  mass  matrix  and  it  is  the  same  matrix  used  in  linear 
structural  dynamics.  Using  the  global  shape  function  defined  in  Eq.  18,  the  mass  matrix  of 
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the  element  can  be  written  as 
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where  m  is  the  mass  of  the  beam  element  and  l  is  its  length.  Note  that  a  consistent  mass 
approach  has  been  used  in  defining  the  mass  matrix.  It  can  be  demonstrated  that  this  mass 
matrix  leads  to  exact  modeling  of  the  rigid  body  inertia,  while  a  lumped  mass  approach  would 
lead  to  a  wrong  modeling  of  the  inertia.  Using  the  expressions  of  the  elastic  and  inertia  forces 
previously  obtained,  in  the  absolute  nodal  coordinate  formulation  the  equations  of  motion 
of  the  finite  element  take  the  following  simple  form: 


Me  +  Qfc  =  Qa, 


(23) 


where  Qk  is  the  vector  of  the  elastic  forces,  and  Qa  is  the  vector  of  applied  nodal  forces. 
While  the  mass  matrix  is  a  constant  matrix,  the  vector  of  elastic  forces  is  highly  nonlinear 
function  of  the  absolute  nodal  coordinates.  The  preceding  equation  can  be  written  as 

Me  =  Q,  (24) 

where  the  vector  Q  =  Qa  —  Q*..  Since  the  mass  matrix  is  constant,  efficient  and  accurate 
numerical  procedures  can  be  used  to  solve  the  preceding  system  of  equations  for  the  vector 
of  the  generalized  accelerations  e.  For  instance,  a  Cholesky  decomposition  of  the  symmetric 
positive  definite  mass  matrix  can  be  made  once  at  the  beginning  of  the  integration  and  used 
throughout  the  entire  numerical  solution. 
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6  PERFORMANCE  IN  DYNAMIC  APPLICATIONS 


Twoproblems  have  been  investigated  in  order  to  compare  between  the  performances  of  the 
corotational  procedure  and  the  absolute  nodal  coordinate  formulation  in  dynamic  applica¬ 
tions.  As  expected,  problems  are  encountered  when  the  corotational  procedure  is  used  to 
solve  some  simple  multibody  applications.  Since  the  ANSYS  code  can  not  be  used  to  model 
complex  multibody  problems,  in  this  section  only  simple  multibody  applications  are  consid¬ 
ered.  These  applications  are:  (1)  free  falling  of  a  pendulum,  and  (2)  non-smooth  motion  of 
a  four  bar  mechanism. 

Pendulum  Problem  The  first  dynamic  problem  considered  in  this  section  is  the  free 
falling  of  a  very  flexible  two  dimensional  beam  under  the  effect  of  gravity.  The  beam  is 
connected  to  the  ground  by  a  pin  joint  at  one  end,  as  shown  in  Fig.  8.  The  beam  has  a 
length  of  1.2  m,  a  circular  cross  section  with  an  area  of  0.0018  m2,  a  second  moment  of  area 
of  1.215E-08  m4,  and  a  modulus  of  elasticity  of  0.700E+06  Pa.  In  the  original  configuration, 
the  beam  is  horizontal  and  has  zero  initial  velocity. 

Two  cases  are  considered  in  the  analysis  of  the  falling  pendulum.  In  the  first  case,  the 
beam  is  assumed  to  fall  under  the  effect  of  gravity,  while  in  the  second  case  the  beam  is 
accelerated  by  increasing  the  gravity  constant  to  50  m/s2.  The  results  of  the  two  models 
of  the  pendulum  are  obtained  using  the  absolute  nodal  coordinate  formulation  and  the 
corotational  procedure  proposed  by  Rankin  and  Brogan  [15].  Three  models  were  considered 
to  simulate  the  motion  of  the  free  falling  pendulum.  These  models  employ  12,  40  and 
100  finite  elements.  The  configurations  of  the  free  falling  pendulum  at  different  time  steps 
predicted  using  the  absolute  nodal  coordinate  formulation  and  the  12-element  model  are 
shown  in  Fig.  9.  Figure  10  shows  the  transverse  deflection  of  the  midpoint  of  the  pendulum 
using  the  three  different  models.  From  the  results  presented  in  this  figure  it  is  clear  that 
the  12-element  solution  leads  to  accurate  results.  The  solutions  obtained  using  the  40  and 
100-element  models  are  identical. 
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In  the  case  of  the  large  value  of  the  acceleration  constant,  the  deformation  of  the  beam 
becomes  much  larger.  The  configurations  of  the  pendulum'at  different  time  steps  predicted 
using  the  absolute  nodal  coordinate  formulation  and  the  40-element  model  are  shown  in  Fig. 
11.  In  this  case,  the- 12  element  solution  does  not  lead  to  very  accurate  results,  due  to  the 
large  deformation.  However,  the  results  obtained  using  40  elements  axe  accurate,  as  shown 
by  Fig.  12  where  the  transverse  deflection  of  the  midpoint  of  the  beam  is  plotted  versus 
time. 

Figure  13  shows  the  transverse  deflection  of  the  mid  point  of  the  pendulum  obtained  using 
ANSYS  when  the  gravity  constant  is  equal  to  9.81  m/s2.  The  12-element  model  does  not 
converge.  Furthermore,  the  40-element  solution  diverges  after  0.7  sec  despite  the  simplicity 
of  the  model.  Only  when  a  large  number  of  elements  is  used,  convergence  is  achieved  using 
the  corotational  formulation.  It  is  important  to  point  out  that  changing  the  number  of  steps 
and  the  number  of  convergence  iterations  does  not  result  in  an  improvement  of  the  results. 
This  convergence  problem  is  attributed  to  the  use  of  local  rotations  as  nodal  coordinates 
in  the  corotational  formulation.  In  this  problem,  the  relative  rotation  between  the  surface 
coordinate  system  and  the  element  coordinate  system  becomes  larger  than  30°  when  a  small 
number  of  elements  is  used.  This  leads  to  problems  when  the  corotational  formulation  is 
used,  as  explained  in  Section  2. 

When  the  gravity  constant  is  increased  to  50  m/s2,  the  corotational  procedure  fails  in 
the  simulation  of  the  motion  of  the  simple  pendulum.  In  this  case,  a  high  value  of  the  ac¬ 
celeration  and  a  relatively  high  mass  produce  large  inertia  forces,  and  this  results  in  large 
deformations  and  large  angular  velocities.  The  high  inertia  forces  and  angular  velocities, 
which  are  characteristics  of  multibody  applications,  pose  serious  problems  when  the  coro¬ 
tational  formulations  are  used.  As  demonstrated  by  the  results  presented  in  Fig.  14,  100 
elements  are  not  enough  to  achieve  convergence,  and  the  solution  diverges  after  0.3  sec. 
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Four  Bar  Mechanism  The  simulation  results  of  the  simple  pendulum  previously  pre¬ 
sented  in  this  section  clearly  demonstrate  some  of  the  serious  problems  that  can  be  encoun¬ 
tered  in  the  simulations  of  very  simple  multibody  systems  as  the  deformation  and  speed 
increase.  In  this  section,  another  multibody  example,  the  four  bar  mechanism  shown  in 
Fig.  15,  is  considered.  The  dimensions  and  material  properties  of  the  links  of  the  four  bar 
mechanism  are  shown  in  Table  3.  All  components  of  the  mechanism  are  made  of  steel  and 
have  a  circular  cross  section  with  diameter  equal  to  0.4  m.  This  system  is  designed  to  obtain 
high  values  of  the  angular  velocities  of  the  connecting  rod  and  the  follower  as  compared  to 
the  angular  velocity  of  the  crankshaft.  In  this  system,  complete  rotations  of  the  crankshaft 
are  possible,  as  the  Grashoff ’s  law  gives: 

s  +  l  =  1.2  <  1.21  =p  +  q,  (25) 

where  s  and  l  are  the  lengths  of  the  shortest  and  longest  links,  and  p  and  q  are  the  lengths 
of  the  other  two  links.  However,  the  difference  between  the  two  sides  of  Eq.  25  is  very 
small,  and  this  makes  the  motion  non-smooth.  In  the  case  of  rigid  body  motion,  the  angular 
velocities  of  the  connecting  rod  and  the  follower  are  presented  in  Fig.  16  as  functions  of 
the  angle  of  rotation  of  the  crankshaft  assuming  a  unit  value  for  the  angular  velocity  of  the 
crankshaft.  It  is  clear  from  the  results  presented  in  this  figure  that  when  the  rotation  of 
the  crankshaft  is  close  to  0,  27r,  47t,  ...  the  angular  velocities  of  the  connecting  rod  and  the 
follower  change  dramatically  in  a  very  short  time.  The  system  is  assumed  to  be  driven  by  a 
moment,  shown  in  Fig.  17  as  a  function  of  time,  applied  to  the  crankshaft,  and  the  effect  of 
the  gravity  force  is  taken  into  consideration. 

Figure  18  shows  the  transverse  deflection  of  the  midpoint  of  the  connecting  rod  predicted 
rasing  the  absolute  nodal  coordinate  formulation.  The  transverse  deflection  is  determined  as 
the  distance  of  the  midpoint  from  a  straight  line  that  connects  the  two  ends  of  the  connecting 
rod.  It  is  clear  from  the  results  presented  in  Fig.  18  that  up  to  approximately  0.75  sec  the 
motion  is  very  smooth  and  the  deformation  of  the  connecting  rod  remains  small.  After  0.75 
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sec,  the  crankshaft  completes  a  full  revolution  and  there  is  a  jump  in  the  angular  velocity 
leading  to  a  large  deformation. 

Several  simulations  have  been  performed  using  the  corotational  formulation  implemented 
in  ANSYS,  using  different  numbers  of  steps  in  the  integration  routine.  In  the  first  simulation, 
20,000  time  steps  were  chosen  while  maintaining  the  option  of  automatic  stepping  active. 
This  simulation  configuration  leads  to  the  solution  shown  in  Fig.  19,  where  the  global  vertical 
position  of  point  A  on  the  crankshaft  is  presented  and  compared  to  the  solution  obtained 
using  the  absolute  nodal  coordinate  formulation.  Before  0.75  sec  there  is  no  difference 
between  the  two  solutions,  but  after  that  the  two  curves  diverge.  It  is  clear  that  for  the 
corotational  procedure  to  converge,  a  smaller  integration  step  is  required.  In  order  to  achieve 
this,  the  automatic  stepping  option  is  removed  in  a  second  simulation.  This  change  improves 
the  results  significantly,  as  demonstrated  by  the  results  shown  in  Fig.  20.  However,  there 
are  still  differences  when  the  deflections  are  considered  instead  of  global  positions  of  nodes, 
as  demonstrated  by  the  results  shown  in  Fig.  21.  Furthermore,  increasing  the  number  of 
time  steps  to  40,000  does  not  lead  to  a  better  improvement  of  the  results,  as  shown  by  the 
results  of  Fig  22. 

In  this  four  bar  mechanism  problem,  the  total  deformation  of  the  bodies  remains  small, 
but  the  angular  velocity  experiences  jumps  each  time  the  crankshaft  completes  a  full  cycle. 
Hence,  in  the  vicinity  of  that  configuration,  the  displacement  increments  are  large  within 
a  single  time  step,  and  the  results  given  by  the  corotational  procedure  are  not  accurate. 
Theoretically,  convergence  can  be  achieved  as  the  time  step  approaches  zero.  This,  however, 
may  lead  to  excessive  error  accumulation  in  practice.  In  fact,  Fig.  22  shows  that  the  con¬ 
necting  rod  has  the  same  pattern  of  vibration  as  previously  predicted  by  the  absolute  nodal 
coordinate  formulation.  Nonetheless,  in  the  case  of  the  corotational  procedure  implemented 
in  ANSYS,  after  about  1  sec,  a  phase  shift  develops,  and  this  shift  cannot  be  corrected  with 
a  further  increase  in  the  number  of  time  steps. 
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7  SUMMARY  AND  CONCLUSIONS 


Many  multibody  applications  are  characterized  by  motion  discontinuities,  high  inertia  forces 
and  high  and  discontinuous  angular  velocities.  In  this  investigation,  two  finite  element  proce¬ 
dures,  the  corotational  technique  and  the  absolute  nodal  coordinate  formulation,  which  can 
be  used  for  the  solution  of  large  deformation  problems,  are  presented  and  their  computational 
performances  are  demonstrated  using  several  numerical  examples.  In  this  investigation,  the 
limitations  of  the  corotational  formulation,  that  has  been  implemented  in  several  general 
purpose  finite  element  codes,  are  demonstrated  when  flexible  multibody  applications  are 
considered.  It  is  shown  that  the  incremental  procedure  can  be  computationally  expensive 
in  large  deflection  problems  as  compared  to  the  non-incremental  absolute  nodal  coordinate 
formulation  in  which  the  nodal  coordinates  are  defined  in  a  fixed  inertial  frame.  The  absolute 
nodal  coordinate  formulation  leads  to  a  constant  mass  matrix  which  is  the  same  as  the  mass 
matrix  used  in  linear  structural  analysis.  Therefore,  the  inertia  forces  are  linear  functions  in 
the  accelerations  and  the  dynamic  equations  of  motion  do  not  include  any  quadratic  velocity 
terms.  The  elastic  forces,  on  the  other  hand,  are  highly  nonlinear  function  of  the  nodal 
coordinates  even  in  the  case  of  linear  elastic  models. 

In  the  case  of  static  analysis  of  beam  large  deflection  problems,  it  is  demonstrated  that 
the  absolute  nodal  coordinate  formulation  leads  to  accurate  results.  On  the  other  hand,  it  is 
shown  that  the  corotational  procedure  can  be  computationally  expensive  and  can  lead  to  a 
lock  in  the  solution  because  of  the  presence  of  the  rotations  in  the  set  of  nodal  coordinates. 
The  performance  of  the  absolute  nodal  coordinate  formulation  in  dynamic  problems  has 
also  been  evaluated  using  a  flexible  pendulum  and  a  flexible  four  bar  mechanism.  Due 
to  the  limitations  on  the  amplitudes  of  the  rotations  in  the  corotational  procedure,  such 
a  formulation  can  fail  in  the  simulation  of  simple  multibody  systems,  as  demonstrated  by 
the  results  presented  in  this  study.  In  applications  characterized  by  high  inertia  forces  and 
motion  and  velocity  discontinuities,  serious  problems  can  be  encountered  in  the  simulation 
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of  flexible  multibody  systems. 

It  was  demonstrated  that  the  results  obtained  using  the  absolute  nodal  coordinate  formu¬ 
lation  agree  well  with  the  results  obtained  using  the  floating  frame  of  reference  formulation 
in  the  case  of  small  deformation  problems  [7,  19].  The  absolute  nodal  coordinate  formula¬ 
tion,  however,  can  be  used  as  the  basis  for  developing  a  new  generation  of  flexible  multibody 
codes  that -can  be  used  in  the  small  and  large  deformation  analysis  of  flexible  multibody 
systems,  as  demonstrated  in  this  investigation. 
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Table  1.  Elastica  problem:  global  rotations  of  the  free  end  node.  Exact  and  ANSYS  solutions 


p/pv 

cr 

1.015 

1.063 

1.152 

1.293 

1.518 

1.884 

2.541 

4.029 

9.116 

Qfree-end  eXaCt 

20° 

40° 

60° 

80° 

100° 

120° 

140° 

160° 

180° 

Q/rec-end  1 0  elements 

33.25° 

41.51° 

59.35° 

79.53° 

99.82° 

119.98° 

130.77° 

- 

Q/rcc-end  20  elements 

33.85° 

42.15° 

59.80° 

79.79° 

99.93° 

119.98° 

131.84° 

- 

- 

4/2 


Table  2.  Elastica  problem:  global  rotations  of  the  free  end  node.  Exact  and  absolute  nodal  coordinate 

formulation  solutions 


PI  P 

cr 

1.015 

1.063 

1.152 

1.293 

1.518 

1.884 

2.541 

gjfgjg 

9.116 

Qfree-end  eXaCt 

20° 

40° 

60° 

80° 

100° 

120° 

140° 

180° 

o  Jree-end  1 0  elements 

21.41° 

38.56° 

78.55° 

98.53° 

118.48° 

138.56° 

175.57° 

Qfree-end  20  elements 

22.53° 

39.74° 

80.04° 

100.17° 

120.20° 

140.22° 

ess 

176.11° 

Table  3.  Parameters  used  in  the  simulation  of  the  four-bar  mechanism 


Body 

m  [kg] 

A  [m2] 

/[m4] 

/[m] 

El Pa] 

Crankshaft 

4.9323 

1.257E-03 

1.257E-07 

0.5 

2.1E+11 

Coupler 

6.9052 

1.257E-03 

1.257E-07 

0.7 

2.1E+11 

Follower 

5.5242 

1.257E-03 

1.257E-08 

0.56 

2.1E+11 

Fig.  1.  Corotational  procedure 


a)  Undeformed  configuration 


a)  Deformed  configuration 


Fig.  2.  Absolute  nodal  coordinate  formulation 
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Fig.  3.  Cantilever  beam  bent  into  a  full  circle  by  an  end  moment.  ANSYS  solution 


Fig.  4.  Rotations  of  the  free  end  node  of  a  cantilever  beam  subject  to  end  moments 
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Fig.  5.  Cantilever  beam  bent  into  a  full  circle  by  an  end  moment.  Absolute  nodal  coordinate 

formulation  solution 


Deflection  (m) 


Fig.  10.  Transverse  deflection  of  the  midpoint  of  the  pendulum  for  different  models.  (a=9.81  m/s2) 


Fig.  12.  Transverse  deflection  of  the  midpoint  of  the  pendulum  for  different  models.  (a=50  m/s2) 


Fig.  13.  Transverse  deflection  of  the  midpoint  of  the  pendulum  for  different  models.  (a=9.81  m/s2) 
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Fig.  14.  Transverse  deflection  of  the  midpoint  of  the  pendulum  for  different  models.  (a=50  m/s2) 


Angular  Velocity  (rad/s) 


0.65  m 


Fig.  15.  The  four  bar  mechanism  in  the  original  configuration 


Angle  of  rotation  of  the  crankshaft  (deg) 


Fig.  16.  Angular  velocities  of  the  connecting  rod  and  the  follower  in  the  case  of  rigid  body  motion 
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Fig.  17.  Moment  applied  to  the  crankshaft 


Fig.  18.  Transverse  deflection  of  the  midpoint  of  the  connecting  rod  obtained  using 
the  absolute  nodal  coordinate  formulation 
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Fig.  19.  Global  vertical  position  of  point  A  on  the  crankshaft 
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Fig.  20.  Global  vertical  position  of  points  on  the  crankshaft 
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